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Using dynamic renormalization group we study the transport in driven diffusive systems in the pres- 
ence of quenched random drift velocity with long-range correlations along the transport direction. 
In dimensions d<4 we find fixed points representing novel universality classes of disorder-dominated 
self-organized criticality, and a continuous phase transition at a critical variance of disorder. Nu- 
merical values of the scaling exponents characterizing the distributions of relaxation clusters are in 
good agreement with the exponents measured in natural river networks. Published in: Phys. 
Rev. E, vol. 58, p. 168 (1998) 
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I. INTRODUCTION 



Various interacting driven systems self-organize into 
critical steady states, optimizing in that way their re- 
sponse as a single functional unit An important 
feature of these systems is their response in the pres- 
ence of disorder. Effects of disorder on critical properties 
in steady states have been investigated ||, jj| in cellu- 
lar automata and coarse-grained continuum models. It 
has been recognized that disorder changes local relax- 
ation rules and breaks some symmetries of the dynam- 
ics, which may result in a qualitatively different global 
dynamic state. A distinct class of physical phenomena 
in driven systems exhibit the scale-free behavior only in 
the presence of disorder. Examples include energy trans- 
port in the integrate-and-fire oscillators with diversity [Q , 
and Barkhausen noise in spatially disordered ferromag- 
nets H . Fluctuations of optimal path in heterogeneous 
materials Q and landscape evolution due to river net- 
works flowing in naturally fractal environment (7|-|^] also 
belong to this class of dynamical systems. 

Most studies of self-organized criticality (SOC) have 
been done on sandpile automata, in which nonlinearity 
responsible for SOC is due to threshold condition of top- 
pling. In the continuum equation of motion for the dy- 
namic variable height h(x,i), this leads to an infinite se- 
ries of relevant operators fi n d 2 h n Recent numerical 
simulations of stochastic automata with " soft" threshold 
reveal new universality classes of SOC and a phase transi- 
tion when probability of toppling is varied |jlT| . Comple- 
mentary to numerical simulations, the renormalization 
group studies of continuum equations are aimed to char- 
acterizing the critical behavior at large distances and long 
times. Hwa and Kardar (HK) [|l2| introduced a trans- 
port equation which is compatible with all symmetries 
of anisotropic flow in open diffusive system, with leading 
nonlinearity \d\\h 2 generated by nonlinear friction. 

In this work we study the transport equation of open 
diffusive systems in the presence of quenched random 
drift velocity. We adopt an anisotropic d-dimensional 



model for the height h(x, t) transport with HK non- 
linearity \d\\h 2 , and introduce quenched disorder via a 
new term p{x)d\\h, which locally breaks the joint inver- 
sion symmetry h — > —h, xu — > — afii. The symmetry is 
globally restored by assuming the distribution of disor- 
der with zero mean, which thus excludes global current 
through the system. We also consider long-range correla- 
tions of disorder along the direction of transport varying 
with distance as ~ jx^ 2+s and weak (anti)correlation in 
perpendicular direction (see below). In our model p(x) 
represents spatially-varying local velocity of profile fluc- 
tuations, which is motivated by mass transport in real- 
istic granular and river flow with a preferred direction 
of drainage. It should be stressed that our model dif- 
fers from continuum models of SOC studied so far both 
in its symmetry properties and in correlations of defects. 
Meanwhile, in models of randomly driven interfaces |l5| ] 
a similar disorder term appears in a physically different 
context. 

Using dynamic RG in the hydrodynamic (HD) limit we 
show that this type of disorder represents a relevant per- 
turbation in dimensions d < 4, leading to a new disorder- 
induced scaling behavior. We calculate the critical expo- 
nents at new fixed points in the e = 4 — d and 5- expan- 
sion to leading order |l4j . It is interesting to note that in 
the absence of disorder the e-expansion to leading order 
yields the exact critical exponents in the HD limit, as 
discussed in detail in Ref. [Q. Using scaling arguments 
eligible for directed dynamic processes which generate 
sclf-affine structures, we also determine the avalanche ex- 
ponents in terms of the anisotropy exponent. 

The organization of the paper is as follows: In Sec. II 
we introduce the stochastic differential equation with dis- 
order and discuss motivation for the long-range disorder 
correlations. In Sec. Ill the details of the dynamic renor- 
malization group analysis are given. In Sec. IV we discuss 
the critical behavior at disorder-induced fixed points for 
various physical values of the parameters and their rel- 
evance for the problems of river networks and strongly 
disordered dynamical systems. A short summary of the 
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results is given in Sec. V. 



II. STOCHASTIC EQUATION OF 
DISORDER-DOMINATED NETWORKS 

We start with the anisotropic diffusion equation for the 
height transport with nonlinear friction 

dh A 

— = u\\djh + v L d]h - -^d\\h 2 -p(x)d\\h + r) , (1) 

and a time-dependent nonconserving Langevin force 

(tj(x, t)r]{x't')) = 2D5 (d \x - x')8{t - t') . (2) 

The new random term proportional to p(x) locally breaks 
joint inversion symmetry, which is obeyed by the remain- 
ing terms . We assume the distribution of p(x) as 



(p(x)) d = ; (p(x)p(x')) d = jf(x - x') 



(3) 



with f{x) — 



'II 



For consistency of the pertur- 



bation expansion (see below) we will chose a = 2 — 5 
and c ~ 0{2 — z). In Eq. ([!]) anisotropy signals the 
existence of a preferred direction of mass flow, which 
is the subject of two nonlinear terms: (\/2)(d\\h 2 ) and 
7/(2;) (<9|| h) 2 . The motivation for long-range disorder cor- 
relations is as follows: We assume that Eqs. (|l])-(||) de- 
scribe the evolution of heights (e.g., of a granular pile 
or landscape) which eventually leads to a self-organized 
structured landscape with a network of channels, along 
which the material is being eroded. It is important to 
keep in mind that these channels appear dynamically as 
a result of diffusion, which is influenced by an interplay 
of the above two nonlinear terms. Therefore, an initial 
configuration that is based mainly on the configuration 
of disorder, helps to imprint the channels by setting lo- 
cally most probable drift paths. However, since the sys- 
tem is open and repeatedly perturbed by the noncon- 
serving noise 77, the once established network of channels 
is likely to evolve under further perturbations, reaching 
a new stationary configuration, in which effects of dis- 
order are altered. An example of dynamically modified 
disorder effects can be found in the field-driven random 
Ising model, in which pinning by local random fields ap- 
pears weakened by sweeping an avalanche of flipped spins 
over certain pinning centers. The size of an avalanche is 
the subject of the dynamics itself. Another interesting 
example is represented by erosion of natural landscapes 
due to water flow. In the course of evolution, the orig- 
inally preferred local drift directions become uneffective 
at sites which are found inside the correlated area that 
already drains to a different direction. Observations in 
natural river basins reveal [ p~5|j3~6| a persistent correlation 
between the average soil slope at a site x and area s(x) 
that drains to that point as X7h(x) ~ [s(x)] _1 / 2 . Here 



the drainage area s(x) is not fixed but is determined self- 
consistently by the dynamics itself. A nice example of 



this relation at work is shown in Ref. |16|, where a pro- 
cedure numerically iterated to self- consistency yields a 
self-similar river network. In the stationary critical state 



<0 



we have s(x) ~ x,, cf>(x±/x'^ 



where D\t is the fractal 



dimension with respect to parallel length, and £ is the 
anisotropy exponent (see below). It is reasonable to as- 
sume that for xu — > 00 the scaling function </>(r) behaves 
as a power of r, i.e., <p{r) ~ r v . Therefore, for the inter- 
mittent dynamic regime [where eroded material from the 
area s{x) is accumulated at point x building up a shear 
stress cr(x), and erupting when the stress exceeds a crit- 
ical value a A , the leading nonlinear term is proportional 



to 



x[. 



{Vh) 2 /x\ 

In order to mimic the above processes in which the ef- 
fects of disorder are being dynamically modified, we only 
fix the disorder correlations in the direction of trans- 
port. The transverse correlations are then determined 
self-consistently by the fluctuations in the critical steady 
state [Q. Notice that the difference 2 — z is a measure 
of the strength of critical fluctuations (see Sec. III). 



III. DYNAMIC RENORMALIZATION GROUP 
ANALYSIS 

The dynamic renormalization group consists of elim- 
inating fast modes with subsequent rescaling: t — > 
b z t, xh — > 6a; 11, x± — » b^x±, h — > b x h, where z, 
C, and x are the dynamic, anisotropy, and roughness 
exponent, respectively. Naive dimensions of the cou- 
pling constants in Eq. (|l|) are then obtained from the 
scaled equation: dh/dt — b z ~ 2 v\\d 2 h + V'^^v^&^h — 
b z+ x -i\ d ^ h 2 + b z - 1 +^pp(x)d ll h + b z -x + ^T], where ac- 
cording to Eqs. we have p, p = — |(a — c£), and 
= — \ [1 + C(d — 1) + z]. In (fc, o;)-space the equation 



(J-T/ — 2 

of motion becomes 

h(k, to) = Go(fc, uj) rj(k,w) 

A f d d q dto' 
!l 2 J (2^p 27 



(4) 



h(q, uj')h{k — q,u — uj') 



d d q dw' . , r _ 

I T^-T — (k\\-q\\)p{q)h{k-q,uj-uj) 



(2ir) d 2ir 



with the propagator Go(k,u>) = v\\k 2 + v^k 2 ^ 



It- 



erating Eq. (|g) and eliminating fast modes leads to a 
diagramatic expansion (see Fig. 1). In the hydrody- 
namic (HD) limit k±_ — > 0, ui — > 0, k\\ <C 1, keeping 
the lowest respective orders of fcii , we calculate the one- 
loop contributions to the recursion relations [£ = log 6, 
S d = 2 1 - d n- d / 2 /T(d/2)}: 



dvu 
~d£ 



z-2 



32 



2irSdW 



(5) 



2 



dD 

dl 



D 



z - 2 X - (d- 1)C - 1 + ^S d w 



d\ 



= A 



\(7-d)-\-^-S d w 



a + c( + nS d w} 



(6) 



(7) 



(8) 



Here the effective couplings u and w are found to be u = 

d-l d—l 

^2 ( — ] and w = 3* [ — ) . A few comments 

are in order: (a) As usual in systems with quenched ran- 
domness, the perturbation expansion is made at fixed 
random noise p(q) , and subsequently the graphs are av- 
eraged over the distribution of p(q) leading to a dashed 
line with a cross, which carries 79||~ a Qj_. All graphs must 
be connected before this step is taken, thus leading at 
most to one dashed-crossed line per loop. Due to the 
quenched nature of the random noise [cf. Eq. (§)], a loop 
with dashed-crossed line involves no frequency integra- 
tion. However, averaging over the dynamic noise accord- 
ing to Eq. (Q) leads to a solid-circled line with two propa- 
gators and a factor 2D, and an integration over the inter- 
nal frequency, (b) According to Eq. (|I|), the wiggly line 
associated with the vertex A/2 carries ik\\, with k being 
the momentum of the incoming line, whereas the wiggly 
line associated with the vertex p(q) carries the momen- 
tum — of the outgoing line. Hence both graphs 
in Figs, la and lb for the renormalized propagator are 
proportional to «(fc|| —q\\)ik\\ and thus do not contribue to 
the vertex i/±. Therefore, we have dv±/d£ = v±[z— 2£], 
leading to £ = z/2. This argument is valid to all or- 
ders in the HD limit, (c) Additional three graphs for A/2 
and 7 (not shown) which are obtained by replacing the 
dashed-crossed line in Figs. Id and If with a solid-circled 
line (there are three such graphs for for A/2 and three 
for 7, corresponding to a circled line along one of the 
three sides of triangle), however, give a null contribution 
(same as in Ref. |l2|]). Similarly, a contribution of the 
graph for dynamic noise D, which is obtained by replac- 
ing the dashed-crossed line in Fig. lc by a solid-circled 
line, vanishes. 

On approaching a fixed point, we have from Eq. (||) 
X = [z(3 -d)-2 + nS d w] /4. From Eqs. (|) -(|) we find 



du 



dw 
~d£ 



9?r „ 9tt 



G4 



S d u 



S d w 



6 - tt-S^u - 3irS d w 
16 



(9) 



(10) 



where the small expansion parameters are e = 4 — d and 
6 = 2— a, and we have chosen c = (2— z)/2 = 1— C- Notice 
that this choice of c is selected by the structure of the true 



expansion parameters u and vo, so that e and 5 appear 
as their anomalous dimensions, respectively. Also, the 
disorder correlations f(xu,x±) become isotropic when 
C = 1, corresponding to the isotropic transport. It should 
be stressed that for anisotropic disorder correlations no 
additional parameters are generated to leading order. 
Eqs. (P)-(p^|) have four fixed points (u*,w*): (G) Gaus- 
sian (0,0), (P) Pure (9^,0), (R) Random (0, 3^ 

and (M) Mixed 



;From Eq. 



the dynamic exponent is obtained as 
„ 3tt 



:-!2 



S d u* - 2nS d w* 



(11) 



Using Eq. (|ll|) and the above scaling relations between 
z, C and X, we find the exponents in the e- and S- ex- 
pansion, which are shown in Table I. Also shown are 
the exponents a and r for the probability distribution 
of duration P(t) ~ t~ a and size of relaxation clusters, 
P(s) ~ s~ r , which can be expressed in terms of ( us- 
ing the following scaling arguments (see also Ref. M). 
For strictly directional diffusion, the clusters can be vi- 
sualized as effectively planar structures with the frac- 
tal dimension D\\ = 1 + (. The average size of clusters 
scales as (s) ~ L de , where L is the linear system size and 
di = 1 for the self-affine clusters (for which C<1)- On 
the other hand, (s) ~ L D ii < - 2 ~ r ', and the scaling relation 
D|| (r — 1) = a — 1 holds in the steady state. Using these 
relations we find r = (1 + 2Q/(1 + £) and a = 1 + £• 



IV. UNIVERSALITY CLASSES OF 
DISORDER-INDUCED CRITIC ALITY 

As seen from Table I, the fixed point (G) represents 
mean-field SOC, which becomes unstable for dimensions 
d<A both with respect to nonlinearity and disorder. Rel- 
ative stability of the other three fixed points in (u, w)- 
plane depends on the parameters e and S and on the 
initial value of the ratio w/u. In Fig. 2 we show the flow 
diagrams of Eqs. C|o|— p~0^) for S = 1, e = 1 and e = 2. (For 
convenience we use reduced couplings U = nS d u/32 and 
W = irS d w.) In the case e = 1, competition between 
disorder and the A-term leads to two different types of 
behavior, which are separated by the line W/U = 1. The 
fixed point (M), whose domain of attraction is the line 
W/U = 1, is unstable in the direction perpendicular to 
the critical line W/U — 1, representing a phase transition 
from pure (HK) to disorder-controlled SOC with increas- 
ing variance of disorder. We find qualitatively the same 
behavior for short-range correlations (S = 2) in d = 2. 
In the case of long-range correlated defects (6 — 1) in 
two dimensions (cf. Fig. 2) the (M) fixed point moves 
to the negative [/-region and becomes spirally attractive. 
The entire first quadrant flows towards the pure HK fixed 
point (P). The flow lines are first attracted to a section 
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of the curve connecting the (R) and (P) fixed points ap- 
proaching fixed point (P) under a nonzero angle. 

Taking the analytic continuation to 8 — ► 1 and e — > 1 or 
2, corresponding to physical d = 2+1 or 1+1 dimensions 
p8[ , respectively, we obtain numerical values of the expo- 
nents, which are listed in Table II. Here 8 = 1 was taken 
as a typical example of long-range correlations. Notice 
that in contrast to e, which is restricted to integer val- 
ues, the parameter 8 may vary continuously in the range 
< 8 < 2. It is noteworthy that the exponents at all 
three fixed points in d = 2+1 are close to values measured 
in natural river networks (RN). For river basins around 
the world the exponents are found 0] as: r = 1.41 — 1.45, 
a = 1.67 - 1.92, C = 0.67 - 0.92, and Hack's exponent 
h = 0.54 — 0.6 satisfying the scaling relation h = 1/a. 
The roughness exponent for large length scales jl9| was 
found in the range x — 0.3—0.55. Variations in the values 
of the exponents depend on geographical location where 
they have been measured, and can be related to locally 
dominated erosion mechanisms 20 1. With regard to the 
results in Table II we would like to point out the follow- 
ing: (i) In the absence of disorder 7 = 0, corresponding 
to the limit studied by Hwa and Kardar in Ref . |L2| , the 
exponents are within the range of the above RN expo- 
nents, indicating that the HK model of flowing granular 
piles captures the basic features of landscape evolution. 
It should be stressed that in this case (7 = 0) the values of 
the exponents are exact, and are not a subject of higher- 
order corrections in the perturbation expansion (see Ref. 
|l2"| for details), (ii) For finite disorder 7 two more 
fixed points are accessible, depending on the initial values 
of 7 and A. Therefore, variations of numerical values of 
the exponents can be attributed to different universality 
classes, which are accessible for varying initial strengths 
of disorder. In particular, a decreases from 1.75 at HK 
limit to 1.72 at (M) and eventually to 1.66 at fixed point 
(R) by increasing disorder 7 at fixed A (see Fig. 2 top). 
In addition, the exponents at fixed points (M) and (R) 
vary with the range of disorder correlations 8 (values of 
8 in the interval < 5 < 2 correspond to long-range 
correlations). For instance, for 8 = 1/3, e = 1 we have 
at (M) fixed point: ( = 0.83, a = 1.83, and r = 1.45. 
It is interesting to note that the same values for cluster 
exponents are obtained by the numerical procedure in 
Ref. |l(|. According to the discussion in Sec. II, we have 
that for £ = 0.83 the exponent in the leading nonlin- 
ear term becomes Du — ((1 — () = 1-69, which is close to 
2 — 8 = 1.67. On the other hand, for shorter disorder cor- 
relations, e.g., for 8 = 3/2, we find a = 1.64 and r = 1.39. 
The importance of disorder for river networks has been 
also pointed out in Ref. plj , where numerical simula- 
tions of a cellular automaton model of randomly pinned 
landscape evolution yields the exponents very close to 
those at fixed point (M) (see Table II). Moreover, in 
the absence of pinning, the same authors [^l[ found the 
exponents close to HK fixed point (P) given in Table II. 



At the (R) fixed point, representing the universality 
class in strong disorder limit (A = 0), the exponents are 
in very good agreement (cf. Table II) with the results of 
numerical simulations of Ref. ||] , suggesting small higher- 
order corrections, (d e = 0.98 + 0.02, £ = 0.66 + 0.02, and 
t = 1.40 ± 0.02) for self-affine networks with disorder- 
dominated basins in two dimensions. It has been ar- 
gued in the literature that the problem of optimal 
path in strongly disordered medium and Eden growth 
processes also belong to this universality class. In these 
systems the disorder effects are dynamically modified. 
Eden growth is not defined as a disordered problem, how- 
ever, an effective quenched disorder with long-range cor- 
relations is self-generated by the blocking effects of pre- 
viously occupied sites (see second citation in Ref. ||). 
Similarly, in the above mentioned example of disorder- 
dominated basins in two dimensions sites that are already 
connected at time t influence the course of the process 
at later time steps. In numerical simulations of cellular 
automata models, such as done in Refs. ||,^l],[ll]], for in- 
stance, a particular range of disorder correlations is not 
specified. The exponents are measured in the emergent 
stationary critical state, which is obtained after many 
successive updates. To our knowledge, a stochastic dif- 
ferential equation for dynamically varying disorder effects 
in these systems has not been considered so far. Here we 
argue that Eq. ([!]) with disorder correlations of the type 

/(S^- 

these dynamical systems. 

Following general scaling arguments of Ref. |lj] we es- 
timate behavior of the order parameter, defined in anal- 
ogy to cellular automata by the average outflow current 

< J(W) >= Jy j d d - l x x j{L h x^t,W). (12) 

In the steady state (J(W)) = 1, and exhibiting fluctua- 
tions near the transition, where we have {J(W)) ~ W 13 . 
For small disorder the local current is j ~ h 2 , thus we 
have j(L h x x ,t,W) ~ tfxjQ)- 1 ^ b^x ±1 b~ z t, b~^W) 



.!•,,, ~Xj_ * might capture the critical properties of 



or j(L h x ±1 t,W) 



2x/C 



4>(tx 



-*/< 



Wx 



-s/C 



) . Insert- 



ing into Eq. (12) and after extracting formally the In- 



dependence we find = [2% + ((d — l)]/8. The directed 
diffusion in our model represents some kind of a contact 
process, therefore for t — > 00 and X|| — > 00 the following 
scaling relation holds: (3/v = a'. Here v is the parallel 
correlation length exponent, and a' = a — 1 is the expo- 
nent of the survival probability distribution. At (M) fixed 
point for 8 — 1 and e = lwe find (3 = 0.78 and v = 1.08. 
Similar value j3 = 0.86 was found in the stochastic cellu- 
lar automaton (see first citation in Ref. |jTT||). 

V. CONCLUSIONS 

We have demonstrated that our transport equa- 
tion with quenched disorder in the drift velocity with 
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anisotropic long-range correlations describes two novel 
universality classes of critical behavior in open diffusive 
systems. For finite disorder we find SOC relevant for the 
scaling properties of fractal river networks. For low dis- 
order a crossover occurs to the asymptotic behavior con- 
trolled by the HK fixed point At critical variance 
of disorder a continuous phase transition occurs between 
the two different types of steady states: channel-type 
flow for strong disorder and low friction, and surface-like 
flow for low disorder and high friction. Comparison of 
the numerical values of the avalanche exponents at fixed 
points in 2+1 dimensions with the exponents measured in 
natural river networks is quite satisfactory. Our anal- 
ysis suggests that natural river networks may result from 
the interplay between quenched disorder and an effective 
nonlinear friction. Variations in the range of disorder 
correlations < S < 2 appear as a possible underly- 
ing mechanism which explains the observed variations in 
the exponents of natural networks. A distinct universal- 
ity class of disorder-induced self-organized criticality is 
represented by the (R) fixed point of our model, where 
A = 0. Evidence collected by numerical simulations in 
Refs. suggest that a number of other disordered dy- 
namical systems should have the same critical behavior 
described by the (R) fixed point. In the present work we 
pointed out the importance of long-range correlations in 
this class of self-organizing disordered systems. 
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TABLE II. Numerical values of the critical exponents at 
various fixed points for S =1, e — 1 (upper part) and e = 2. 
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FIG. 1. One-loop diagrams with nonzero contributions to: 
(a,b) renormalized propagator, (c) dynamic noise D, and ver- 
tices (d,e) A/2, and (f) 7. The symbols are defined in the 
bottom line. 
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FIG. 2. Flow diagrams for 5 = 1 and (top) e = 1 and (bot- 
tom) e = 2. Large circles represent fixed points described in 
the text. 
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